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Abstract. Trigonometric polynomials are widely used for the approximation of a smooth func- 
tion from a set of nonuniformly spaced samples. If the samples are perturbed by noise, a good 
choice for the polynomial degree of the trigonometric approximation becomes an essential issue to 
avoid overfitting and underfitting of the data. Standard methods for trigonometric least squares 
approximation assume that the degree for the approximating polynomial is known a priori, which is 
usually not the case in practice. We derive a multi-level algorithm that recursively adapts to the least 
squares solution of suitable degree. We analyze under which conditions this multi-approach yields 
the optimal solution. The proposed algorithm computes the solution in at most 0{rM + Af^) oper- 
ations (A/ being the polynomial degree of the approximation and r being the number of samples) by 
solving a family of nested Toeplitz systems. It is shown how the presented method can be extended 
to multivariate trigonometric approximation. We demonstrate the performance of the algorithm by 
applying it in echocardiography to the recovery of the boundary of the Left Ventricle of the heart. 
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1. Introduction. The necessity of recovering a function from a finite set of 
nonuniformly spaced measurements arises in areas as diverse as digital signal pro- 
cessing, geophysics, spectroscopy or medical imaging. The measurements {sj}'j^i are 
often distorted by several kinds of error. Hence a complete reconstruction of the 
function from the perturbed data s| = sj + i/j is not possible. Often the function 
to be reconstructed is smooth, in which case a trigonometric polynomial of relatively 
low degree (compared to the possibly huge number of samples) can provide a good 
approximation to the function. This trigonometric approximation may be found by 
solving the least squares problem 



mm 

P&Pm 



^\pix,)~s'^\'w,, (1.1) 



where Wj > are weights and Pm is the space of trigonometric polynomials of degree 
less than or equal to M. 



Many efficient algorithms have been developed to solve (1.1), e.g., see the articles 



p2| , 1^, ^ |ri|, 10 1 . But surprisingly little attention has been paid to the problem of 
how to control the smoothness of the approximation in order to avoid overfitting and 
underfitting of the data. An adaptation of the smoothness of the approximation can 
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be achieved of instance by providing a suitable upper bound for the degree M of the 
space Pm in (1.1). In most of the aforementioned algorithms a necessary requirement 
to get useful results in applications is that a good a priori guess of the degree of the 
trigonometric approximation is available. However a priori it is not clear what is a 
suitable degree for the solution, in terms of how to choose a reasonable degree M 
when solving (^]^) . Determining M by "trial and error" is certainly not a satisfactory 
alternative. 



It is the goal of this paper to derive an efficient algorithm that computes the 
trigonometric approximation which provides the "optimal" balance between fitting 
the given data and preserving smoothness of the solution. Here optimality is meant 
in the sense that the solution has minimal degree among all trigonometric polynomials 
that satisfy a certain least squares criterion. The algorithm recursively adapts to the 
least squares approximation of optimal degree by solving a family of nested Toeplitz 
systems in at most 0{Mr + M^) operations, 



If the data {sj}^^]^ were (i) unperturbed and (ii) stem from sampling a trigono- 



metric polynomial (with degree less than r/2), then the solution of (1.1) would auto- 
matically have the appropriate degree, since the original function could be completely 
recovered in this case. However the assumptions (i) and (ii) are rarely met in ap- 
plications and controlling the smoothness of the solution becomes essential to avoid 
overfitting and underfitting of the data. If wc choose the upper bound for the degree 



in (1.1) too large, the solution will almost always take on the maximal possible degree, 



hence being too wiggly and picking up too much noise (overfit), see also Figure LI 
(a)-(b). In the extreme case 2M + 1 = r we will get an interpolating polynomial, 
mostly with strong oscillations and far away from approximating the function between 
the given samples. On the other hand, if we choose AI too small, then the approxi- 



mation will be very smooth but poorly fitting the given data (underfit). Figure 1.1(c) 
illustrates this behavior. The "regularized" trigonometric approximation obtained by 
the algorithm proposed in this paper - to which we will refer as Levins on- Galerkin 



algorithm - is shown in Figure 1.1(d) 



The paper is organized as follows. In Sections § and ^ we present the main 
results, including the Levinson-Galerkin algorithm and a theoretical analysis that 
clarifies under which conditions this algorithm provides optimal results. In Section ^ 
we show how properly chosen weights can be used as simple but efficient tool to 
precondition the least squares problem. Some aspects of extending the algorithm to 
multivariate trigonometric polynomials are discussed in |[ In Section ^ we present 
some applications in echocardiography. 



Before we proceed we introduce some notation and conventions. The inner prod- 
uct is denoted by (•, •), and the conjugate transpose of a matrix A by A* . The space 
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(a) Original signal and perturbed samples (b) Least squares approximation using a 

too large upper bound for the degree of 
the polynomial (overfit) 
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(c) Least squares approximation using a 
too small upper bound for the degree of 
the polynomial (underfit) 



(d) Regularized approximation by pro- 
posed Levinson-Galerkin algorithm 



Fig. 1.1. Controlling the smoothness of the solution is essential for trigonometric approxima- 
tion from perturbed data in order to avoid overfitting and underfitting of the data. The proposed 
Levinson-Galerkin algorithm automatically adapts to the least squares solution of optimal degree. 



of trigonometric polynomial of degree equal to or less than M is defined as 



= lp:p{x)= J2 Cfee^"'''^ . 



(1.2) 



The norm of p{x) = J2h=-M CfeC^'^''^^ e Pm is given by 

1 \ ^ / M 



= 1 1 \p{x)\''dx\ =(^J2 l^^l') = 



(1.3) 
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where c — {ck}^L_^j. In some applications it is advantageous to deal with complex- 
valued polynomials (see also Section hence we do not restrict ourselves to the case 
of real-valued trigonometric approximation. 

For a — [a_j\,/, a_M+i- • • ■ , ciM-i, o-m] G C^*^"''^ we define the orthogonal projec- 
tions P/v by 

Pata = [0, . . . , 0, a_Ar, a^N+i: • • • , aN-i,aN,0, . . . , 0] (1.4) 
for N = 1,2, .... M and identify the image of Pjv with the 2N + 1-dimensional space 

Let p^^ and be trigonometric polynomials of degree M and N respectively, 
with coefficients vectors c^^ e C2^'^+\c^ € 0^^+^ If < M, then we can always 
interpret p^ as polynomial of degree M by adding to appropriate number of zero- 
coefRcients and by doing so we are embedding the vector into a zero-padded vector 
of length 2M + 1. Wc will henceforth tacitly assume that such an embedding has been 
made, when we compute expressions such as ||c*^ — c^ll- 

2. Multi-level least squares approximation. A standard method in numer- 
ical analysis to find the optimal balance between fitting the given data and preserving 
smoothness of the solution is to introduce a regularization parameter. The best value 
of this regularization parameter is then determined for instance by generalized cross 
validation jist or via the L-curve [ pO[ . Here wc understand regularization not as a way 
to stabilize ill-conditioned problems, but in a broader context as a means of finding 
the best compromise between fitting a given set of data and preserving smoothness 
of the solution. As we will see, in our case it is not necessary to introduce an addi- 
tional parameter, since we can regularize the smoothness of the solution by varying 
the parameter M of the space Pm in which we are searching for the solution of (^]^) . 

For the derivation of the algorithm we consider first the following situation. As- 
sume p4x) = E^)'=-A'.(c,)fce2"'==' G Pn, and let s" = {Sj^};^i,2A/ + 1 < r with 
Sj — Sj + Vj = {xj ) + Vj be given noisy samples satisfying 

U^sr<e\\s^r. (2.1) 

For convenience we assume that r, the number of samples, is odd. 

The aim is to approximate from the data {Sj}^^]^. Let us first assume that we 
already know that we are searching for our least squares solution in the space Pn, ■ In 



this case the coefficient vector of the polynomial that solves (1.1) is the least squares 
solution of 

WVm,c = Ws'', (2.2) 
where Vjv, is a r x (2iV* + 1) Vandermonde matrix with entries 

{VM),,k = e'"*'-"^- , ]^l,...,r,k = ~N^,...,N, (2.3) 
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and W = dia,g{{ y/wj}). 

We will discuss the role and specific choice of the weights in more detail in Sec- 
tion ^. To reduce the notational burden we absorb the weight matrix W in the 
Vandermonde matrix and in the sampling values. Thus for given degree, M say, we 
consider the linear system of equations 

Va/c = (2.4) 

where Vm is now the r x (2M + 1) "weighted Vandermonde" matrix. We will denote 
the least squares solution of ( ^.4D by c^^^' = {c\!'^^}kL~M ^^^'^ the corresponding 
polynomial is p(^^^\x) = Efcl-M ci^'^^e^-*'^-. 

Since in general we do not know the optimal degree or level M of the space Pm in 
which we should solve the least squares problem, the situation becomes considerably 
more complicated. If we want to solve (1.1) under the information (2.1) without 
knowing the degree of the polynomial, one may argue that we have to accept any 
trigonometric polynomial p(x) = X^fcL-w "^fc^^'^*'^^ with ||V;vc — s^H < e\\s^\\ as an 
approximate solution to p,, since it is compatible with the only knowledge we have 
on the data. 

In general there may be infinitely many such polynomials, which raises the ques- 
tions of how to find a polynomial p that yields a small approximation error \\p^ — p\\ 
and at the same time can be computed efficiently. 

2.1. A multi-level algorithm and an efficient stopping criterion. The 

heuristic considerations above suggest the following approach. 

Algorithm 1. Set N = and solve Vqc^°^ ~ s^. //c^") satisfies the condition 
IjVoc^'^-' — s^W < e||s'^||, take c'*'-' as solution. Otherwise set N = N + 1 and solve 

Vnc^""^ = , (2.5) 

until c^^^ satisfies for the first time the stopping criterion 

li^jvcW-s^ll <£||s1|. (2.6) 

at some level N = Nq. Set c'^^°^ = c'^^\ The approximation to is then p^^°'>{x) = 

(No) 2TTikx 

The stopping criterion (|2.6|) is well-defined, since it is definitely satisfied for N = 
(r— l)/2, in which case the left side in (2.6) equals 0. Thus Algorithm |l| selects among 
all least squares solutions p*-^-*, iV = 0, . . . , (r — l)/2 that polynomial with minimal 
degree. 

Algorithm |l| and stopping criterion ( |2.6| ) can be justified by the following theo- 
retical considerations. 

One readily verifies that the matrices Vn, N = Q, . . . , (r— 1)/2 satisfy the relations: 
(i) there exists a left-inverse such that 

V+Vn = In, with V+ = {V^VNy^V^, (2.7) 

where In is the identity matrix on C^^+^. 
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(ii) Let a G C^*^+^ be the coefficient vector of some p E Pm- Then 

VNa = VMa for aU > M, a e 0^*^+^ (2.8) 

In (ii) we have made use of the fact that the coefficient vector a can be interpreted 
as coefficient vector of a polynomial of degree N by extending it to a vector of length 
2N + 1 via zero-padding. The matrix-vector multiplication Vmo, and equation ( |2.8| ) 
should be understood in this sense. 

Lemma 2.1. If N > thenc^^'> satisfies |1 -s'^]] < £||s^||, hence stopping 

criterion (2.6) always becomes active at some level Nq < N^. 

Proof. Note that V^V^ is the orthogonal projection into range(VAr) and s e 
range(VArJ ^ range(VAr) for < N, hence V^V^s = s. Therefore 

II VivcW _ =\\VnV+{s + i.) - (s + = - VNV+,yf (2.9) 
= M\'-\\VNV+,yf <e\s^\\^ (2.10) 

where we have used condition ( ^.l| ) in the last step. It follows from ( 2.1Cl| ) that 
Algorithm |l| terminates at some level Nq < N^.. □ 

The following lemma shows that from the viewpoint of numerical stability it is 
advisable to keep the level N of the space Pj\i in which we search for our solution as 
small as possible. 

Lemma 2.2. cond{V^VN) > cond{VljVM) for N > M. 

Proof. Since 

Pm{V^Vn)Pm = V^iVm for M < N, 

Cauchy's Interlace Theorem |l6| implies that cond(V^VAr) > cond{V^,iVM) for N > 
M. □ 

In the sequel we demonstrate that the fact that Algorithm terminates at some 
level < A'^* is a desired property in many cases. We show that stopping criterion ( |2.6| ) 
is even optimum in a number of cases. 

Let us first consider two special cases: (i) noisefree samples and (ii) uniformly 
spaced samples. 

2.1.1. Noisefree samples. Any reasonable stopping criterion has to satisfy the 
following lemma. 

Lemma 2.3. For noisefree data the stopping criterion (2.6) yields the exact so- 
lution. 

Proof. One readily verifies that Algorithm |l| terminates at level . Hence for 
N = N,: 

IM^^ ~p4 = -c4 = WVj^s^ -c4 = Wv^v^c^ - c.|| = o, (2.11) 

since V^^Vnc^ — c^ for N > N^. U 



Lemmas 2.3 and 2.2 together show that stopping criterion (2.6) yields the opti- 
mum solution for noisefree data while providing maximum numerical stability. 
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2.1.2. Uniformly spaced samples. If the sampling points Xj, j — 1, . . . ,r are 

uniformly spaced and we choose wj — 1/r as weights then a simple calculation shows 
that Vat is unitary on C^^+\ i.e., V^Vn = /at for = 0, 1, . . . , (r - l)/2. 
In this case 

||c, - cW|| = ||c, - V^s'W = ||c, - Vr^VN,c, - V^iy\\. (2.12) 
N > implies V^Vn, = In and hence 

l|c*-cW|| = llF>||. (2.13) 

Note that 



(2.14) 



since VatV^ is an orthogonal projection. Equation (2.14) yields 



for M < . 



Consequently 



c^*''^ II < II c*-cW II ifiV, <Af<7V. 



(2.15) 



(2.16) 



Thus for uniformly spaced samples any stopping criterion should terminate Algo- 
rithm 1^ at the latest at iV = iV* . Under a mild condition on the coefhcients c* we 
can show that the proposed stopping criterion provides the optimal solution among 
all least squares solutions. 

Proposition 2.4. Assume that the samples are regularly spaced. Then the solu- 
tion computed via Algorithm |^ satisfies 



lb* -P^^°'|| < b* -P^^^ll for all N>N, 
If furthermore satisfies 

\\{In, - Pn)c4 > \\{In. - PnW^^i^W 

then 

Ib^-p'^^^'ll < lb*-p''^'|| for all N. 



(2.17) 



(2.18) 



(2.19) 



Condition ( 2.18| ) is satisfied e.g., if all coefficients of are larger than the relative 
noise level, i.e., |cfc| > e||s^||- 

Proof Lemma pj| yields that Nq < N^, thus ( ^.ie| ) implies ( |2.17| ). 
To prove assertion (^.19 ) we only have to show that 

lie - c(^«)|| < ||c, - c(^)|| for aU N < 
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For N < note that (c* — V^Vat.c*) is orthogonal to V^iy, since 
={VN,c,,V;^V^iy) = {VN,c,,VNVj^i^), 

hence 



(2.20) 
(2.21) 



{c,-v;^Vn,c,,v^v) = 0. 

Therefore 

lie, - = lie - V+Vn,c, + = lie, - V+Vn,c4' + \\VM^ . 

In order to prove |lc, — c'^°'|l < |lc, — c'^^^H for aU N < N^, we need to verify 
lie, - V+Vn,c4^ + \\V^iy\\^ > \\V;^^J^f. Since 



\\c.-V+Vn,c4^ ^ J2 l(c,)fcP - |l(/jv-Pjv)c,|p 

lA;l=Af+l 



and 



\k\=N+l 

the result follows now from the assumption ( ^.1^ ). □ 



Remark: Proposition |2.4| shows that the least squares polynomial that gives the 
best approximation to is not necessarily of degree iV, . 

2.1.3. Noisy nonuniform samples. For noisy nonuniformly spaced data we 
observe that 



lb. 



-P^^)|| = lk* 



-cW| 



^^^iv.c,|| + ||y+Hl, 



and for iV > iV, 



(2.22) 



IIP* ~P' 

since ||c, — V^Va^.c, |1 = in this case. 

If Vat is not unitary then |1 V^;/|l is not necessarily monotonically increasing with 
increasing level N. One can argue heuristically that since ||V^|| is increasing with 
increasing level N due to Lemma |2.2| , we may fairly assume that ||V^i^|| will also 
increase (although not strictly monotonically). Also from the viewpoint of numerical 



stability it is reasonable to keep the level N small, since by Lemma 2.2 we know 



that cond(V^VAr) > cond(y^jVA/) for N > M. This together with (|2.22|) suggests to 
choose a stopping criterion which terminates at or before level TV, , which is guaranteed 



for stopping criterion (2.6) by Lemma 2.1 



We can conclude that the stopping criterion will provide excellent results if the 
noise level e is small or if the condition number of V^Vn is small (which implies that 
Vn is approximately unitary) . In order to verify the latter it is useful to have estimates 



for the condition number of V^Vn- We will address this issue in Proposition 3.1 
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2.2. A Toeplitz system and trigonometric approximation. Instead of di- 
rectly solving Va/c''*^-' = it is more efficient in our case to consider the normal 
equations 



(2.23) 



The reason is that from a numerical point of view the structural properties of the 
matrix VIjVm are much more attractive than those of Vmi which in turn leads to 
faster numerical algorithms, see also Section ^. 

Set Tm — V^jVm then a simple calculation shows that the entries of the hcrmitian 
matrix Tm are 



M)k,l 



^2'Ki{k — l)xj 



kj 



(2.24) 



Tm is a Toeplitz matrix, since the entries {TM)k,i depend only on the difference k — l. 
Obviously Tm is invertible if 2M + 1 < r. 



Following result is just a reformulation of (2.23) together with relation (2.S), but 
since it plays a key role in Section ^ it is helpful to state it in detail (cf. also ||l^ ) . 

Theorem 2.5. Given the sampling points Q < xi < . . . ^Xr < ^, samples {s^jj^x; 
positive weights {wj}^^^^ and the degree M with 2M + 1 <r. The polynomial p^'^'^^ G 



Pm that solves (1.1) is given by 



p^^'\x) 



M 

E 

k=-M 



^(M) 2Trikx p 

-■■m ^ ^ -*■ J 



M ■ 



where its coefficients c[,*^^ satisfy 



Tmc 



(2.25) 



(2.26) 



with 



l(M) 



27Tikxi 



for \k\ < M, 



and Tm as defined in (2.24) 



(2.27) 



3. Weights as simple preconditioner. Vandermonde matrices are known to 
be ill-conditioned, if the nodes xj are clustered To improve the stability of the 



systems (2.4) and (2.26) we can use the weights as simple diagonal preconditioner. 
This leads to the problem of how to choose the weights wj . 

We propose to use the size of the area of the Voronoi region |Q associated with 
the sampling point Xj as weight Wj . In 1-D this reduces to 
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This choice is motivated by the foUowing observations. 

In this section we let Vjv denote the Vandermonde matrix defined in (2.3) without 
weights. Let p £ Pm with coefficient vector c. Since 

r 

(Tmc, c) - {WVmc, WVmc) = WWVMcf - ^ \p{xj)?Wj, (3.2) 
the inequahty 

r 

Cillcf <^b(x,)pu;, <C2|lcf (3.3) 

holds for all c G ^sM+i constants Ci = Ami„ and C2 = Amaa;, where Xmin and 
Amaa; dcnotc the minimal and maximal eigenvalue, respectively, of Tm- 

(i) The lower bound of Tm is mainly determined by the large gaps in the sam- 
pling set. Suppose there is a large gap in the sampling set and denote the interval 
corresponding to this gap by T (hence Xj ^ F). Choose a trigonometric polynomial 
p £ Pm which, like the prolate spheroidal functions, concentrates most of its energy 
in the interval T . Then the sampling values of p will not pick up any information 
about the main concentration of the polynomial energy. Consequently if we use no 
weights (or set Wj — 1) we get 

Y,\p{x,)\''^\\pf = \\cf. 



For such a sampling set the lower frame bound Ci in the inequality (3.2) must be 
small. Generically, large gaps and the ensuing lack of information always results in 
bad condition numbers. This problem cannot be fixed by preconditioning. 

(ii) On the other hand, we can choose a trigonometric polynomial that is mainly 
concentrated in the region where the sampling points are located. In this case the 
same local information is counted and added several times. Thus 



b(^.)P»IWP = l|c||' 



and the upper constant C2 in ( |3.3| ) will be large. Yet, as mentioned in (i) a cluster 
will not contribute much to the lower bound and to the uniqueness of the problem. In 
this case the condition number is large, because too much local information is given 
in certain areas of the polynomial. 

Problem (ii) can be addressed by introducing properly chosen weights. The idea 
is to compensate for the local variation of the sampling density by using weights in 
inequality ( |3.3| ). Suppose that < xi < X2 < • • ■ < 2;^ < 1 is a sampling set in [0, 1]. 
Then a natural choice for the weights is Wj — [xj^i — Xj^i)/2. Thus if many samples 
are clustered near a point Xj , then the weight Wj is small. If Xj is the only sampling 
point in a large neighborhood, then the corresponding weight is large. This choice 
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has not only been confirmed by extensive numerical experiments , but also by the 
following optimization approach. 

A standard approach for the construction of prcconditioners for a matrix A is the 
following. One attempts to find the matrix P in a given class A4 of matrices (e.g., 
the class of all circulant matrices or the class of all diagonal matrices) which solves 

mm\\I - PA\\f , (3.4) 
PeM 

where ||.||f denotes the Frobenius norm. 

In our setting this translates to the following optimization problem 

mmlll ~ {WV)*WV\\f , (3.5) 

where V is the class of all r x r diagonal matrices and / is the {2M + 1) x (2M + 1) 
identity matrix. 



Note that we require that Wj > whereas (3.5) could in principle yield weights 
that violate this condition. However since we will make use of ( p.5[ ) in our actual 
algorithm, we are somewhat sloppy here. 

An alternative approach is to consider the solution of 

inin{cond[{WV)*WV]} 
subject toW eV. 

This optimization problem can be transformed to a general eigenvalue problem, see |^ , 
which can be solved by convex optimization algorithms. 

In the simple case of regular sampling it is easy to check that the solution of both 
optimization problems is given by W = diagdy/wj}) with Wj = (xj+i — Xj^i)/2 = 
1/r. However in the more interesting case of nonuniform sampling neither prob- 



lem (3.5) nor (3.6) does in general have an analytic solution. Thus using these ap- 
proaches for the actual construction of a preconditioner would be ridiculous, since 
the computational costs to solve these optimization problems are considerably larger 



than solving the trigonometric approximation problem. Nevertheless, solving (3.5) 



and (3.6) numerically for a variety of different examples is useful to get insight in the 
type of weights obtained by these approaches. 

The numerical results confirm the choice of the Voronoi-type weights defined 



in (3.1). Sampling points in densely sampled areas are assigned a small weight, 
whereas sampling points in sparse sampled regions are assigned a large weight. Two 
typical comparisons of the weights obtained via optimization and the Voronoi weights 
are illustrated in Figure ^. In the first case we consider a sampling set with high 
density at the endpoints and strongly decreasing density towards the center. The 
weights obtained by solving ( |3.5| ) and ( |3.6D are almost identical and are very close to 
the Voronoi weights, as can be seen in Figure ||(a). The difference at the endpoints is 
probably due to boundary effects. 
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Fig. 3.1. Comparison of weights obtained by different approaches. 

In the second example we consider a random sampling set with several areas 
with high sampling density and relatively few samples between these clusters. Again 
all three approaches give weights that show a similar behavior, see Figure §(b). The 
condition number of the non- weighted Toeplitz matrix in this example is 33, compared 
to the significantly smaller condition number 3.3 when using Voronoi-type weights. 



Using the weights obtained via (3^) gives cond(rA/) = 3.1, and for the weights 
resulting from ( |3.6| ) we get cond(TA/) = 2.9, which is only a slight improvement 
compared to the Voronoi-type weights. 

Obtaining good estimates for the condition number of a Toeplitz matrix is a 



difficult problem. It is gratifying that by using the weights defined in (3.1) it is 
possible to get an upper bound for the condition number. 

Proposition 3.1 (Grochenig, |]l8|). Assume that the sampling set {xj}^^^ sat- 
isfies 

max(2;j+i - Xj) — 7 < ^ 
and set Wj ~ (xj+i — Xj-i)/2. Then the condition number of the Toeplitz matrix Tm 



defined in ( |2.24| ) is hounded by 



niTM)<[l^y. (3.7) 



4. A Levinson-Galerkin algorithm for trigonometric approximation. 

The method described in Algorithm |l| can be seen as a Galerkin-type approach, 
since we try to determine an approximation by searching for a solution in a finite- 
dimensional space spanned by orthogonal polynomials, and by increasing the dimen- 
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sion of the space we increase the resolution of our approximation by adding more and 
more details. 

When we use Levinson's algorithm to solve ( 2.26 ) for M = 0, 1, . . . , TVo the 
total computational effort would be of 0{Nq), since the solution of each system 
Tmc'^^^ = b'^^^ requires 0{M'^) operations. Using one of the fast Toeplitz algo- 
rithms H, H reduces this effort to 0{kM\ogM) for each level M, where k is the 
number of iterations, thus leading to a total of 0(fciVQ log A^g) operations. In this 
section we show that the systems Tmc*^*^^ = h'^'^'^\ M = 0, 1, . . . , A^o can be solved in 
0{Nq) operations and the total effort (including the calculation of the entries of Tm 
and the evaluation of the stopping criterion ( |2.6| )) for computing p'-^o) jg 0{rNo+ Nq) 
operations. 

Following observation is crucial for the derivation of the proposed Levinson- 
Galerkin algorithm. 

Lemma 4.1. For fixed degree M andM + l letTM,b^^''^ and TM+i,b^^'+'^'> be the 
Toeplitz matrices and right hand sides as defined in ( ^.24 ) and (2.27), respectively. 
Then Tm and b'^^'^^ are embedded in Tm+i and in the following way: 





to 




*2(M+1) 




b-{M+i) 


Tm+1 — 




Tm 






tiM) 




t2(M+l) 




to 




bM+1 



(4.1) 



Proof. (4.1) follows immediately from the definition of Tm and fe*-*^-* and (2.S). □ 
Unfortunately the solutions c^*^) and 0'*^+^' of the systems Tmc^^^^ = b^^'^^ and 
T'm+ic^*^^^-' = are not related is such a simple manner. But we can exploit 

the nested structure of the family {Tm}^,j^i by solving the systems Tmc^^'^^ = 
recursively via a modified Levinson algorithm. The standard Levinson algorithm can- 
not be applied directly, since it only addresses Toeplitz systems, where the principal 
leading sub-matrix and the principal leading sub-vector of the right hand side stay 
unchanged during the recursion, which is not the case here. For Tm+i it does not 
matter, if we enlarge Tm by appending new entries below or above, whereas the right 
hand side fe*^*^-* cannot be rearranged in such a way, the principal leading subvector 
of the right hand side will be changed if we switch from 6*^^^^ at level M to 6'^^^+^^ at 
level M + 1. 

To adapt Levinson's algorithm to our situation, we have to split up the change 



from the system Tmc'-^^ = b'^^"'^ at level M to the system Tm+ic 



(M+l) 



6(^+1) at 



level M + l into two separate steps. Instead of indexing the matrix Tm and the vectors 
^(*^) cf*^) by the degree M, it is therefore advantageous to index them according to 
their dimension. For clarity of presentation we reserve the subscript (M) for the 
degree of the polynomial and its coefficient vector respectively, and use the subscript 
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{£) when we refer to the dimension of the corresponding coefficient vector in C^. Thus 



for even 



[b-i 



and for odd £ we set ' = \b_e-i . . be-i ] 

i i 2 2 

(whence fe'"'^-' ~ bo), analogously for cS^K Further it is useful in the sequel to denote 
t*^^^ = [ti, . . . , tiY" . Then the Toeplitz matrix of size £ x £ is generated by the vector 
[to, [t^^-^'^Y]^ with tk = WjC^''^''''^ according to ( |2^ ). 

Assume we have already solved the system Tmc^'^^^ — fc*-*^-' at level M (with 
£ = 2M + 1) and now we want to switch to the next level M + 1. As we have agreed, 
we do this in two steps. In the first step {£ £ + 1) the Toeplitz system can be 
written as 



{t^'^fEe to 
where Ei is the rotated identity matrix on 

'o 

















bi+1 




2 - 




2 - 



(4.2) 



Ef, - 



1 



I.e., 







System (4.2) can be solved recursively by the standard Levinson algorithm ||2l|, [16| . 
To be more detailed, assume that we have already solved the system T^c*^^-' = 6^^' for 
£ — 2M + 1 and assume further that the solution of the -^-th order Yule- Walker system 
Tiy'^^^ = — 1(^) is available. Then the solution of ( [1.2| ) can be computed recursively by 

Vi+i ^ ibi±i - [t^^YEiC^^^)/Pt 



where 



/3, = to + [t^'-'^VW^ = (1 - a,_iS^)/3,_i 



,(^+1) 



Now we can proceed to the second step 
the Toeplitz system can be expressed as 



to 


(<(^+l))*' 




V g+l 




b_i±i 




Te+1 











I ^ £ + 2 = 2{M + 1) + 1), where 



(4.3) 



with c(^+^) ~ , (t;(^+^))^]^ — rS^'^^^\ Observe that (p~^) cannot be transformed 



to a system of the form (4^) by simple permutations, i.e. just by interchanging rows 



REGULARIZED TRIGONOMETRIC APPROXIMATION 



15 



and columns. Since we have already solved the systems T^+ic'-^^^-' = 6(^+1) and 
T^+iy(^+i) = -^(^+1) we can write 



and 



where we have used in the last step that Ti = [Tg]* which implies that 
is real and therefore to + = to + = 



Note that at each level M we have to check if the stopping criterion (2.6) is 
satisfied. The evaluation of the expression 

r 

Y,\p^'^\x,)^s^^\'w, (4.4) 

can be considerably simplified and by avoiding the evaluation of at the nonuni- 
formly spaced points xj we can reduce the computational effort from 0{Mr) to 0{M) 
operations. 

To do this we define the subspace TZ = {{p{xj)}^^i : p € Pa/} C C with the 
weighted inner product (y, z)-ji = Vj^j'^j UiZ G C. The solution of the least 

squares problem (1.1) is the orthogonal projection of the vector {s^jj^^ G C onto TZ 



and therefore must satisfy 



({p(^^)(x,)} - {.5}, {p^''Hx,)}}-,^ = Y.(P^''H^j) - s^^)p(M)(:c,)w, = 
which implies 

{{p^''Hx,)},{s',})n - {{p^''H^j)},{P^''\^s)}U- (4.5) 

Since 



5: 14 -pW(a:,)pz., =E 141'-. - 2Re(.^{p(*^)(a;,)})^ + b(^^)C 
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by (4^), and because 



M X / M 



j=l j=l ^ m=-M ' ^ n=-M 

M M 



m=~M n=-M ^ j = l 

(Tmc(*-^\c(*^)) = (fe(^^\c(*^)), (4.6) 



it follows that 



E \s^,-p^''Kx,rw, = E 141'-. - {b^''\c^''') ■ (4.7) 

Since X]j=i k^P^j to be computed only once at the beginning of the algorithm, 
the evaluation of (4.4) can be carried out in 0{AI) operations. 

Summing up we have arrived at the following algorithm to compute p^'^'^K 

Algorithm 2 (Levinson-Galerkin algorithm for trigonometric polynomials) . Let 
the sampling points {xj}^^^, sampling values {sj}^^i, weights Wj > and the data 
error estimate e be given. Then the trigonometric polynomial determined in 

Algorithm ^ can be computed in OItNq + Nq) operations by the following algorithm. 



Initialize: to = J2'j=i'^3^ti = Z]j=i ^^'j^^''"' > = J2'j=iSj'Wji<^ = ^"j-^ Is^A^w, 



2/(1) = -ti/to, c(i) = bo/to, /3o = to, ao = -ti/to, ei = {a - bl/to)/a, 1=1. 



vifhile ei > e 

/?£ = (1 - ai^iai-.i)[3i-i 
if £ = 1 mod 2 



Vl+l 



f3e 



2 _ 

bi±i 

2 - 

elseif ^ = mod 2 
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V_ e 



6_| - (cW,tW) 
~ We 



+1) 



_i = |a-(fe(^+i),c(^+i))|/a 



end 



Pe 



y 



e+i) _ 



te+i 



end 

No - i/2 



No 



P 



(JV( 



Remark: Usually one evaluates the final approximation on regularly spaced grid 
points, hence the last step of the algorithm can be realized by a Fast Fourier transform. 
The most costly steps are the computation of the entries of t'-^^ and b'-^\ According 
to Corollary 1 in the entries of Tm and 6*^*^-' can also be computed via FFT 
by embedding the Xj into a regular grid (since the Xj can be stored only in finite 
precision). In this case one automatically gets all entries to,. . . ,tr at once. However 
this trick is only useful if the number of points of the regular grid is of the same 
magnitude as the number of sampling points. Alternatively one may use the numerical 
attractive formulas of Rokhlin or Beylkin for a fast evaluation of trigonometric 
sums at unequally spaced nodes. 

Algorithm || can be simplified for real- valued data, this modification is left to the 
reader. 

Fast Vandermonde solvers require 0{Mr) operations for the solution of Va/ c^*^'' = 
s'^ , cf. psf. It is not clear however if these algorithms can utilize the nested structure 
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of the sequence of matrices {Vm}m in order to give rise to an efficient implementation 
of Algorithm |l|. Moreover it is an open problem if the Vandermonde solvers can be 
extended to multivariate trigonometric approximation. We will see in the next section 
that the extension of Algorithm ^ to higher dimensions is straightforward. 

5. Multivariate trigonometric approximation. An advantage of the pro- 
posed approach, besides its numerical efficiency, is the fact that it can be easily ex- 
tended to multivariate trigonometric approximation. In this section we briefly discuss 
some results for the two-dimensional case. 

We define the space of 2-D trigonometric polynomials P^j by 



M 



M 

3,k=-M 



2-Ki{jx-\-ky) 



(5.1) 



To reduce the notational burden, we have assumed in (5.1) that p has equal degree 
M in each coordinate, the extension to polynomials with different degree in each 
coordinate is straightforward. 

For an arbitrary sampling set {{xj,yj)Yj=i & [0,1)^ and given degree M the 



system matrix according to the 2-D version of Theorem 2.5 is 



(Tm) 



k,l 



kA = 0, 



. 2M . 



(5.2) 



One can easily verify that Tm is a hermitian block Toeplitz matrix with 2M + 1 
different Toeplitz blocks of size (2M -I- 1) x 2M + 1, cf. |2^. For a given sampling 
set let Tm be the block Toeplitz matrix for degree M and Tm+i the block Toeplitz 
matrix for degree M + 1. There is a similar relationship between Tm and Tm+i 
as in the 1-D case. More precisely, denote the Toeplitz blocks of Tm and Tm+i by 
{BM)k,k = 0,...,2M, and {BM+i)k,k — 0,...,2M + 2, respectively. Then one 
readily verifies the following embedding: 



'Af+l = 



{Bm+i)o 




{Bm+i)2{m+i) 
to 



B 



M+1 



{Bm+i)*2(m+1) 




{Bm+i)g 



(5.3) 



to 



(5.4) 
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In ^ Levinson's algorithm has been extended to general block Toeplitz systems. 



With this extension and relation (5.3) at hand, we can easily generalize Algorithm Q 
to 2-D (and along the same lines to multivariate) trigonometric approximation. 

The analysis of the stopping criterion (j^^ in Section ^ can be applied line by 
line to the 2-D (actually to the n-D) setting. The only difficulty arises in the search 
for simple criteria for the invertibility of the block Toeplitz matrix Tm ■ The condition 

(2A/ + l f <r 

is necessary in dimension d > 1, but no longer sufficient, since the fundamental 
theorem of algebra does not hold in dimensions larger than one. In Grochenig 
has derived estimates for the condition number of Tm in higher dimensions. In 2-D 
these estimates can be stated as follows. 

Let Ds{a,b) = {{x,y) G : (x - a)^ + {y - < S"^} be the disc of radius i5 
centered at (a, b). We say that a set {{xj,yj),j = 1, . . . , r} is 5-dense in [0, 1] x [0, 1], 
if {J^j^iDs{xj,yj) 3 [0,1] X [0,1]. In other words, the distance of a given sample 
{xj, yj) to its nearest neighbor {xk, yk) , k ^ j is at most 26. 

Analogously to Section || we choose the size of the Voronoi region Vj associated 
with Xj as weight Wj in the computation of the block Toeplitz matrix Tm in (5.2). 



Suppose that the samphng set {{xj,yj),j = l,...,r}C [0, 1] x [0, 1] is (5-dense and 

^<P^- (5-5) 
Grochenig |l^ has shown that under these conditions 

4 

(2 - ^M5)2 

In particular, for arbitrary 5-dense sampling sets, the block Toeplitz matrix Tm is 
invertible and the 2-D version of Algorithm ^ is applicable. 

5.1. Line- type nonuniform sampling in 2-D. In the following we consider 
a special case of trigonometric approximation in two dimensions. This case arises 
when a function is irregularly sampled along lines. A typical example is illustrated 



cond(TM) < ,^ (5.6) 



in Figure 5.1. Such sampling patterns are encountered for instance in geophysics and 



medical imaging, see also Section 3.2 



Corollary 5.1. Let p e f and let {xj, yj.k}, j — 1, . . . , r, fc = 1, . . . , be a 
sampling set in [0, 1)^ such that 

AiWpW^ <Y^\p{y^^,)\^ <Bi\\p\\^ Ai,Bi>0 (5.7) 

k=l 

for every p G Pm and for all j . Further assume that {xj}j£z is a sampling set such 
that 

r 

A2\\pr < J2 \Pi^^)\' ^ B^Wpf- A2,B2 > (5.8) 
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X 



Fig. 5.1. Line-type nonuniform sampling set 



for every p G Pm ■ Then 

3 = 1 fc=l 

for every p G P^j . 

If {xj} and {yj^.k} o,i"e sampling sets with supi.{xk+i — Xk) = 82 < 
supj,fc(yj,fc+i -Vjm) = 5i < ^^en A; = (1 - 6if,Bi = (1 + 5if,l = 1, 
condition number of the block Toeplitz matrix Tm is bounded by 



k{Tm) < 



{1 + 6,)^ (1 + 62)' 

(l-5l)2(l-<52)2 



Proof. Let x be fixed. Then y ~^ p{x, y) G Pm and hence for all j 

1 J.. 1 

M / \p{xj,y)\'^dy<^\p{xj,yj^u)?<Bi / \p{xj,y)\'^ dy 
^ I. — 1 



by assumption (5^). It follows that 



r ^ r ^ 

i=i J=i 

Now let y be fixed. Then x p{x, y) G Pm and 

1 r 1 

M j \p{x,y)\^ dx <Y,\pi^3^y)? < B2 j \p{x,y)\^dx. 
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Since 



r 1 1 r 

[ \p{x,,y)\'dy^ fY,\p{xj,y)\^dy, (5.14) 



assertion (p.9[) follows by combining (5.12) and (5.13) with (5.14). The estimate of 



the constants Ai , Bi and of the condition number of the block Toeplitz matrix 



follow from Theorem 2.5. □ 



The proof of Corollary dA is due to Grochenig [|T^ . Corollary does not only 
guarantee that p G P^.j can be recovered from its samples p{xj, Vj.k), it provides more. 
An immediate consequence is, that it can be reconstructed by an efficient algorithm 
relying on an successive application of Algorithm and the Gohberg-Semencul rep- 



resentation of the inverse of a Toeplitz matrix. See Section 3.2 for more details and 
an application in medical imaging. 

6. Curve and surface approximation by trigonometric polynomials. 

Trigonometric polynomials can be used to model the boundary or the surface of 
smooth objects. Let us consider a two-dimensional object, obtained e.g. by a planar 
cross-section from a 3-D object and assume that the boundary of this 2-D object is a 
closed curve in M . We denote this curve by f and parameterize it by /*(w) — (^u, Vu)^ 
where Xu and y„ are the coordinates of / at "time" u in the x- and y-direction respec- 
tively. Obviously we can interpret / as a one-dimensional continuous, complex, and 
periodic function, where x„ represents the real part and y„ represents the imaginary 
part of f{u). It follows from the Theorem of Weierstrass (and from the Theorem 
of Stone- Weierstrass |2^ for higher dimensions) that a continuous periodic function 
can be approximated uniformly by trigonometric polynomials. If / is smooth, we can 
fairly assume that trigonometric polynomials of low degree provide an approximation 
of sufficient precision. 

Assume that we know only some arbitrary, perturbed points Sj — {xuj i Uuj ) — 
f{uj) + 6j, j — 1, . . . ,r oi f , and we want to recover / from these points. By a slight 
abuse of notation we interpret Sj as complex number and write 

Sj^Xj+iyj. (6.1) 

We relate the curve parameter u to the boundary points Sj by computing the distance 
between two successive points Sj_i, Sj via 

Ul = (6.2) 
Uj = Uj-i + dj (6.3) 

dj = ^ [xj - Xj-if + [yj - yj-if (6.4) 

for j — 2, . . . , r. Via the normalization tj = tj/L with L = Ur -\- djq we force all 
sampling points to be in [0, 1). Other choices for dj in ( |0| ) can be found in |^ in 
conjunction with curve approximation using splines. 
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Having carried out the transformations ( |6.l| )-( |6^ , we can solve the problem of 
recovering the curve / from its perturbed points Sj by Algorithm |^. 

6.1. Object boundary recovery in Echocardiography. Trigonometric poly- 
nomials are certainly not suitable to model the shape of arbitrary objects. However 
they are often useful in cases where an underlying (stationary) physical process im- 
plies smoothness conditions of the object. Typical examples arise in medical imaging, 
for instance in clinical cardiac studies, where the evaluation of cardiac function using 
parameters of left ventricular contractibility is an important constituent of an echocar- 
diographic examination These parameters are derived using boundary tracing 

of endocardial borders of the Left Ventricle (LV). The extraction of the boundary of 
the LV comprises two steps, once the ultrasound image of a cross section of the LV 



is given, see Figure 6T(a)-(d). First an edge detection is applied to the ultrasound 
image to detect the boundary of the LV, cf. Figure |6.1| (c). However this procedure 
may be hampered by the presence of interfering biological structures (such as papillar 
muscles), the unevenness of boundary contrast, and various kinds of noise |2^. Thus 
edge detection often provides only a set of nonuniformly spaced, perturbed boundary 
points rather than a connected boundary. Therefore a second step is required, to 



recover the original boundary from the detected edge points, cf. Figure 6.1(d). Since 
the shape of the Left Ventricle is definitely smooth, trigonometric polynomials are 
particularly well suited to model its boundary. 



After having transformed the detected boundary points as described in ( |6.lD -(|6.4|) 
we can use Algorithm ^ to recover the boundary. The noise level S depends on the 
technical equipment under use, it can be determined from experimental experience. 
Figure |6.2| (a)-(b) demonstrate the importance of determining a proper degree for the 



approximating polynomial. The approximation displayed in Figure 6.2(a) has been 



computed by solving (1.1) where M has been chosen too small, we obviously have 
underfitted the data. The overfitted approximation obtained by solving ( |l . l|) using 
a too large M is shown in Figure ^^(b) . The approximation shown in Figure |6.1| (d) 
has been computed by Algorithm ^, it provides the optimal balance between fitting 
the data and smoothness of the solution. 

6.2. Boundary recovery from a sequence of images. In cardiac clinical 
studies one is more interested in the behavior of the Left Ventricle over a period of 
time rather than in a single "snapshot" . Thus for a fixed cross section we are given 
a sequence of ultrasound images (usually regularly spaced in time) describing the 
variation of the shape of the LV with time. One cycle from diastole (the state of 
maximal contraction of the LV), passing systole (the state of maximal expansion) to 
the next diastole consists typically of about 30 image frames. Since the behavior of 
the LV is (at least for a short period of time) almost periodic, one can model the 
varying shape of a fixed cross section of the LV as distorted two-dimensional torus, 
which in turn can be interpreted as 2-D trigonometric polynomial. Clearly we have 
to use a different degree for the time coordinate r and for the spatial coordinate u. 
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(a) 2-D echocardiography (b) Cross section of Left Ventricle 




(c) Detected boundary points (d) Recovered boundary of LV computed 

by Algorithm 

Fig. 6.1. The recovery of the boundary of the Left Ventricle from 2-D ultrasound images is a 
basic step in echocardiography to extract relevant parameters of cardiac function. 

Due to interfering biological structures and other distortions it sometimes hap- 
pens that some of the image frames cannot be used to extract any reliable boundary 
information. Thus we have to approximate these missing boundaries from the infor- 
mation of the other image frames. To be more precise, assume that an echocardio- 
graphic examination provides a sequence of ultrasound images Ir taken at time points 
r — 0, 1,...,T — 1, where T is approximately the length of one diastolic cycle (the 
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Fig. 6.2. The approximation in the left image results from using a too small polynomial degree^ 
the approximation in the right image from a too large degree for the trigonometric approximation. 



time points could also be nonuniformly spaced). Assume that some of the images Ir 
provide no useful information, so that we can only detect boundary points {sj^kYk^i 
from the images It-^. , where {TjYj=i is a subset of 0, 1, . . . T — 1. In order to get a 
complete description of the LV for the time interval [0, T], we have not only to ap- 
proximate the boundaries fj from each Ij , but we also have to recover the boundaries 
corresponding to the missing images. In other words we look for a 2-D trigonometric 
polynomial € of appropriate degree AI that satisfies p{Tj,Uj^k) = ixj,k,yj,k) 
where the parameter u is related to Sj,k — Xj,k + iyj.k by formulas (6^)-(6^). This 
approximation can be computed by the 2-D version of Algorithm ^ as indicated in 
the beginning of Section ^. 

Under certain conditions we can use the 1-D version of Algorithm || instead of 
its 2-D version. As long as the assumptions of Corollary |5.l| are satisfied, we can 
compute S P]Ij by a successive application of Algorithm ^ We first approximate 
the boundaries fj for each j separately from its samples {sj.felfc^^^i which yields j 
different polynomials e Pm^- Having done this, the next step is to recover 

the missing boundaries at those time points where no information is available. We 
proceed by approximating successively the missing information "line by line" . We 
choose w = 0, say, and approximate the missing information from the samples ) (u) 
taken at the time points Tj,j — 1, . . . , r. 

Note that the Toeplitz matrices of the systems {Tm)uC^u ^^ — b^u ^^ coincide for all 
u, since the sampling geometry is constant along the w-coordinate (because we have 
recovered all samples at each Tj). Thus we have to solve multiple Toeplitz systems 
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with the same system matrix but different right hand side, ft is weli-fcnown that 
this can be done efficientiy by exploiting the Gohberg-Semencul representation of the 
inverse of the Toephtz matrix |fj]. fn our context this reads as follows. We solve 

(TM)„cr)=e) (6.5) 

for one u by Algorithm^. We can solve now all other systems efficiently by establishing 
{Tm)^^ in the Gohberg-Semencul form 

(r^,)-i = - C/W[C/W]*) /z„ (6.6) 

where L'-^^^ is a lower triangular Toeplitz matrix with z = [zq, zi, . . . , Z2mV i^s 
first column, is an upper triangular Toeplitz matrix with [0, zi, . . . , Z2m]'^ as its 

last column, z being the first column of (Tm)~^- The matrix vector multiplications 
to compute = {TM)u^b^^ can now be carried out quickly using the Fast Fourier 
transform by embedding L^*^' and [/'^^^^ into circulant matrices. 

7. Miscellaneous remarks. For sampling sets with large gaps it can happen 
that the system Ta/c*^^^^ — gets ill-conditioned with increasing degree M and 
therefore Algorithm || may become unstable In this case one can use a different, 
more robust approach, which however comes at higher computational costs We 
solve the system Tmc^*^' = 6^*'^^ iteratively, e.g. by the conjugate gradient method 
until a certain stopping criterion is satisfied at iteration fc, say, yielding the solution 
c^*^-*. We use this solution as initial guess at the next level M + I by setting Cq^^^"'^-' 
[0 (cjj^^^')"^ 0]^. The crucial point in this procedure is to find a stopping criterion that 
guarantees convergence of the iterates, see for more details. 

The computation of the entries of the Toeplitz matrix in Section |^ involves the 
nodes Uj which in this particular case depend on the (perturbed) samples Sj. There- 
fore not only the right hand side but also Ta/ is subject to perturbations. Hence 
in principle one might use the concept of total least squares (see ||T2|) instead of a 
least squares approach. A detailed discussion of this modification is beyond the scope 
of this paper. 
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